Excitations of a Bose-condensed gas in anisotropic traps 
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We investigate the zero-temperature collective excitations of a Bose-condensed atomic gas in 
anisotropic parabolic traps. The condensate density is determined by solving the Gross-Pitaevskii 
(GP) equation using a spherical harmonic expansion. The GP eigenfunctions are then used to solve 
the Bogoliubov equations to obtain the collective excitation frequencies and mode densities. The 
frequencies of the various modes, classified by their parity and the axial angular momentum quantum 
number, m, are mapped out as a function of the axial anisotropy. Specific emphasis is placed upon 
the evolution of these modes from the modes in the limit of an isotropic trap. 



I. INTRODUCTION 



The observation Jl|-|j| of the collective modes of the Bose condensate in ultracold trapped atomic gases has stimulated 
a number of calculations of the collective excitations in these systems. Most of these calculations have been performed 
' using the standard Bogoliubov equations for T — 0-II1J], assuming all the atoms to be in the condensate, although 
, some finite-temperature calculations [jl^,[l3|, which make use of the Hartree-Fock-Bogoliubov equations within the 
t-H ' Popov approximation, have also appeared. Both approaches have been used to investigate the excitations in model 
isotropic traps, and in anisotropic traps typically corresponding to the experimental trap of the JILA group 0. 

At the level of the Bogoliubov approximation, the collective excitations are determined equivalently by solv- 
ing the coupled Bogoliubov equations, the linearized time-dependent Gross-Pitaevskii (GP) equation or a pair of 
hydrodynamic-like equations for the condensate density and velocity field. The methods_of solution have included 
analytical solutions within the Thomas-Fermi approximation for the condensate I 
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15(1, variational solutions of 



the time-dependent GP equation gj9[] and expansion techniques for the solution of the Bogoliubov equations using 
harmonic oscillator bases pJTofl. In this paper we develop an alternative method of solution for arbitrary anisotropic 
traps which is based on the construction of the GP equation eigenstates in terms of a spherical harmonic expansion. 
The expansion of the Bogoliubov quasiparticle amplitudes in terms of these functions then leads to a simplified eigen- 
value problem which is used to map the collective excitation spectrum throughout much of the anisotropy parameter 
space. In addition, by analyzing the mode densities, we are able to provide a more detailed discussion than previously 
available of the evolution of the modes from the isotropic limit. 



II. THEORY 

The ground state properties of a trapped atomic Bose gas are well-represented by the stationary GP equation 

n 2 w 2 



+ V ext (r) +gn c (r) 



$o(r)=^o(r), (1) 



2m 

where n c (r) = |(f>o( r )| 2 is the condensate density normalized to the total number of particles N. V ex t(r) is the external 
confining potential and g = AirTi 2 a/m is the interaction strength determined by the s-wave scattering length a. The 
ground state eigenvalue /i is identified with the chemical potential of the condensate. In the following, we assume 
that V ext (r) is axially symmetric, V ext {r) — V ext (r,6), in which case <J> (r) — &o{r,0). 

Our approach to the solution of the Bogoliubov equations is a basis-set expansion method which makes use of the 
eigenfunctions of the ground state GP Hamiltonian. These solutions can be chosen to be eigenfunctions of L z with an 
angular dependence of e lm ^ ', where m is an integer. To construct these solutions, we convert (Q) to a matrix problem 
by making use of a set of normalized basis functions VVtim(r) = Rnii^Yimi^) which are the eigenfunctions of the 

Hamiltonian Hq — — — r- Vo(r), where Vo(r) is some spherically symmetric potential. One possible choice for this 
potential is the spherical average of the external potential which we assume to be of the harmonic form 

V ext (r, 9) = ~mw 2 (x 2 + y 2 ) + \mio 2 a z 2 . (2) 
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Here, uj r and u> a are the radial and axial harmonic frequencies, respectively. We now expand the external potential 
in terms of Legendre polynomials 

V ext {r 1 9)=Y J V^ t {r)P l { C os9) : (3) 
i 

where 

V%t{r) = ^1 J d9 8inm(cos0)V ext (r,e). (4) 

The I = component is the spherical average V^\(r) = ^mr 2 H) 2 , where u) 2 = ^{2lo 2 + uj 2 ) is the arithmetic mean of 
the squares of the axial and radial frequencies. The only other term in the expansion of the external potential is the 
I = 2 component V^ t {r) = ^mr 2 f3uj 2 where 

2cj 2 a - 2qy 2 

P = ^fT^ (5) 

defines the anisotropy parameter. It varies over the range —1 < f3 < 2, where f3 = —1 corresponds to the infinitely 
long cigar-shaped trap and f3 = 2 the infinitesimally thin pancake-shaped trap. The JILA trap |ij has (5 — 1.40 while 
the MIT trap @ has (3 = -0.991. 

Alternatively, we could use the I = component of the total effective potential V(r) = V ex t(r) + gn c (r) appearing 
in (Q). Regardless of the choice, we define V^r) = Vo(r) + AV(r) where AV(r) is the nonspherical perturbation. 
(Depending on the choice of Vo(r), A7(r) may in fact include an I — component.) Expanding an arbitrary solution 
of (|lj) as 0(r) - E 

nim a nim'4 ; nim(^), the expansion coefficients are determined by the matrix equation 

{£ni - s)a n im + (nlm\AV\n'l'm')a n rii m > = 0, (6) 

n'l'm' 

where e is a possible eigenvalue and e n i are the eigenvalues of ho- For an axially symmetric potential, 

(nlm\AV\n'l'm') = (nlm\AV\n'l'm)5 mm > (7) 

so that states with different m-values remain uncoupled. However the nonspherical perturbation does have the effect 
of coupling basis states with different ^-values. An explicit expression for the potential matrix element is 

/>oo 

(nlm\AV\n'l'm) = VA™, / dr r 2 R n i{r)R n , v {r)AV T {r) (8) 
I Jo 

where AVj(r) are the angular components of AV(r) defined in analogy with those of V ex t(v). The numerical coefficient 
in (|8|) is 



A av = \l^^-(U'00\l0)(U'0m\lm) (9) 



where {hh'mi'mi^z'mz) is the usual Clebsch-Gordon coefficient [Q. If Ay has reflection symmetry in the x-y plane 
as assumed, only states of the same parity are coupled (I and I' both even or both odd). This restriction implies that 
the different m-states can be chosen to have a well-defined parity n = (— 1)'+" 1 with respect to reflections in the x-y 
plane. 

The condensate wave function, <i>o(r) = y/N(f)o(r), is determined by the lowest energy even-parity solution of (^|) in 
the m = subspace and is given by 

$oM) - VN^a^Rru^Yioir) . (10) 

nl 

We can now evaluate the condensate density n c (r, 6), which has an expansion similar to that of the potential: 

n c {r,6) =J2ni(r)Pi (cos 6). (11) 
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Comparing (O) to the square of (10), we obtain 



n i( r ) = £ E y/( 21 + W' + mi'OOmi^a^Rn^r)^^ . (12) 

nl 
n'l' 

These radial functions provide what is needed to complete the specification of the nonspherical potential Ay. Since 
Eq.(^|) depends on the condensate density through AV, this equation must be iterated until a self-consistent solution 
for the condensate wave function is generated. 

We determine the collective excitations of the condensate using the Bogoliubov equations, which are equivalent to 
solving the linearized time-dependent GP equation. As shown in Ref. |T^ |, the Bogoliubov equations can be cast into 
the form 

h 2 ^\r) + 2ghn c (r)ijj ( r ] '(r) = Ef^~\r) 

h 2 4 + \r) + 2gn c (r)h4 + \r) = E?^ + \r) , (13) 

where the functions are related to the quasiparticle amplitudes by ^-^{v) = w;(r) ± Wj(r). The Hamiltonian h 
appearing in these equations is the ground state condensate Hamiltonian shifted by the ground state eigenvalue /i. 
Since these equations are uncoupled, either can be used to determine the excitation energies E^. 

To solve (HI) for tp^ + (r) we proceed as in Ref. jL2| and introduce the expansion ip^ + \r) = J2 a c«Va(r) where the 
4> a (r) are the eigenfunctions of h determined from (^|) , and the expansion coefficients c« are required to be normalized 
as ^f^SaCa^ Ca = EiSij . Substituting this expansion into the equation for (r), we obtain 

+ £ «<W> epcf = Efc$ . (14) 



Within a particular m-subspace, the matrix M a p is given by 



M af} = 2g J drcf>* a (r)n c (r)ct> (r) 

= E ai ni a^l(nlrn\2gn c \n'l'm) (15) 



iln'V 



where off is the a-th eigenvector of Eq. (^|) . Since the matrix M a p is diagonal in both the azimuthal quantum number 

m and the parity II, these quantum numbers also serve to classify the collective modes. Once the eigenvector c« has 
been determined, the density of the i-th mode is given by 



6m(T) ex Mr)4~\r) = ^ (r)^-fc^ tt (r) . (16) 



III. RESULTS 



As an illustration of the technique we consider the case of 2000 rubidium atoms contained within an axially 
symmetric harmonic trap of varying anisotropy. For the JILA trap, uo r /2-K — 75 Hz and uo a /2-K = 212 Hz, giving an 
anisotropy parameter of (3 = 1.40, and an averaged harmonic frequency of u)/2ir — 137 Hz. We keep the latter fixed in 
our calculations and vary the parameter j3. To complete the parameter specification, we take m( 87 Rb) = 1.44 x 10" 25 
kg for the mass of the atoms, and an s-wave scattering length of a ~ llOao = 5.82 x 10~ 9 m. Throughout, lengths 
and energies are expressed in terms of the characteristic oscillator length of the isotropic trap d — {Ti/mu)) 1 / 2 = 
9.21 x 10~ 7 m and the characteristic trap energy HQ = 9.03 x 10~ 32 J, respectively. When expressed in these units, 
the dimensionless condensate wave function depends only on the dimensionless parameters f3 and 7 = Na/d. 

In Fig. 1 we show the excitation spectrum obtained from the solution of (14]). In these calculations the basis 
functions, generated numerically, are limited in number to between 100 and 200, the actual number being controlled 
by the value of the high-energy cutoff set in the basis function expansion. This number provided sufficient accuracy 
over the range of the anisotropy parameter shown. However as /3 increases towards its extreme values of —1 and 2, 
increasingly more basis functions are required to obtain accurate results and our calculations are necessarily curtailed. 
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The modes shown here are those corresponding to axial quantum numbers m = — 4 for both even and odd parity. 
The modes of even parity are shown in Fig. 1(a) and those of odd parity in Fig. 1(b). Let us first examine the 
modes of the isotropic trap along the line (3 = studied previously 0,0 . In order of increasing frequency, we find 
a doubly degenerate mode corresponding to the I = 1 excitations, followed by the triply degenerate I — 2 mode, 
the quadruply degenerate I = 3 mode and a nondegenerate I = mode. (In this description, we do not include the 
additional degeneracy associated with the sign of m which is not lifted by axial anisotropy.) It is on the evolution of 
these modes that we concentrate in the following. 

The lowest doubly degenerate I = 1 modes for the spherical trap correspond to the centre of mass modes of the 
harmonic potential JTfl]. At /? = both modes have a frequency of exactly 1 in units of u>. Of these two modes one 
is the odd parity, m = mode, while the other is the even parity m = 1 mode. (Recall that parity refers to the 
reflection symmetry in the x-y plane.) Both of these modes correspond to a rigid oscillation of the condensate density 
in the axial and transverse directions, respectively. As (3 begins to deviate from zero, the degeneracy of these two 
modes is lifted, one oscillating at the frequency uj a (odd parity) and the other at co r (even parity). When expressed 
in units of u, these modes disperse according to Lo a /ui = + (3 and u r /ui = yl — (3/2. We have plotted these exact 
results for the center of mass modes in Fig. 1 to illustrate the accuracy of our numerical calculations. Any deviation 
of the numerical results from the exact values would indicate the need to increase the number of basis functions in 
the expansions. 

For the experiments in the regime where the axial confining potential is stronger than in the radial direction (such 
as for the JILA trap), the collective modes of interest, other than the centre of mass modes, are those originating 
from the I — 2 mode in the isotropic limit. In particular, the modes observed experimentally jjj are the m = and 
m = 2 even parity modes originating from the I — 2 mode, the lower frequency mode being the m = 2. This latter 
mode has quadrupolar character, with a density fluctuation which is concentrated in the x-y plane. As a result, the 
oscillation is mainly influenced by the radial trap frequency to r and the mode frequency decreases monotonically with 
increasing (3, similar to the behaviour of the m = 1 centre of mass mode. In fact, as will be explained in more detail 
shortly, these particular m — 1 and m — 2 modes are the first two members of a family of modes identified by the 
number '1' in Fig. 1(a). 

We now examine the even parity m = mode originating from I — 2 in more detail. The mode density defined by 
( p"6| ) is shown in Figs. 2 and 3 for a variety of (3 values. Since the mode density is of the form Sn(r) oc f(r,6), the 
interesting spatial dependence is revealed by considering a plane (e.g. the x-z plane) containing the axis of the trap. 
Fig. 2 gives the behaviour of the density along the x and z directions while Fig. 3 gives a contour representation. In 
the isotropic limit, this mode corresponds to the quadrupole I = 2, m = mode for which an expansion (contraction) 
in the radial direction is accompanied by a contraction (expansion) in the axial direction, as shown in the (3 = panel 
of Fig. 2. The contour representation in Fig. 3(b) shows nodal lines making angles of ±54.7° with respect to the 
z-axis. This representation is particularly useful since it indicates the direction of particle flow. In the Thomas-Fermi 

approximation j|, the velocity field is given by = — gVSn which implies that the direction of particle flow is 

normal to the density contours. Although not exact, this relationship between velocity and mode density is still a 
good guide in the Bogoliubov approximation. Fig. 3(b) thus indicates that for this mode there is a circumferential 
flow from the equatorial x-y plane to the polar regions. 

However, as the magnitude of (3 is increased from zero, the character of the mode changes from quadrupolar to 
one which is more accurately described as 'breathing- like'. For negative (3, the cloud breathes in the axial direction, 
while for positive (3, it breathes in the radial direction. Thus, in Fig. 3(a) for (3 — —0.5, we see nodal surfaces 
which are approximately planes perpendicular to the z-axis, indicating a flow of atoms in the axial direction. This 
mode is analogous to the lowest standing wave resonance in an open-ended organ pipe of length L, with wavelength 
A ~ L. Conversely, for (3 — 1.4 in Fig. 3(c), we see a cylindrical nodal surface and a flow which is predominantly 
radial. In view of this behaviour, the mode frequency would be expected to depend mainly upon u> a for negative 
(3 and upon uj r for positive (3, resulting in a dispersion to zero in the (3 — ► — 1 and (3 — ► 2 limits. This is precisely 
the behaviour exhibited by this mode which is the lowest mode labeled '3' in Fig. 1(a). Similar arguments can be 
used to understand the dispersion of the higher frequency modes once a contour representation of the mode density 
is available. For example, we show in Figs. 4(a-c) the contour representation of the lowest odd-parity mode labeled 
'4' in Fig. 1(b), which originates from the I = 3, m = mode in the isotropic limit. Fig. 4(a) for (3 = —0.5 is clearly 
the next standing-wave resonance following the fundamental mode illustrated in Fig. 3(a). Likewise, Fig. 4(c) for 
(3 = 1.4 illustrates a hybrid mode which involves both the radial breathing motion of Fig. 3(c) and an oscillation in 
the axial direction. As a result of this axial motion, this mode has a finite limiting frequency for (3 — > 2. 

Also of interest within the anisotropy regime of the JILA experiment is the m = 3 mode with even parity, originating 
from the I — 3 mode for the isotropic case (the third curve from the bottom, labeled '1' in Fig. 1(a)). The excitation 
energy of this mode is only slightly greater than that of the experimentally observed m = mode and like the latter, 
disperses to zero for (3 — > 2. However, as indicated by the common label of T' in Fig. 1(a), this mode is the third 
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member of the family mentioned previously, which is distinct from the family labeled '3' that the to = mode belongs 
to. The observability of this to — 3 mode would of course depend on the possibility of inducing experimentally an 
excitation having a e 3 ^ azimuthal dependence. 

Although the frequency spectrum in Fig. 1 as a function (3 is quite complex, it is apparent that certain patterns 
have emerged. These correspond to the families labeled by 1 through 4 in the figure, three of which were alluded to 
previously. These families correspond to the modes very recently calculated analytically in the Thomas-Fermi limit 
by Ohberg et al. juj and independently by Fliesser et al. (Some of these mode frequencies were obtained earlier 
by Stringari |Qj.) Their analytic expressions for the mode frequencies of the four families, converted to our notation, 
with Ui corresponding to the i-th group as labeled in Fig. 1, are given below: 

uo'l = (l-/3/2)mw 2 

lo 2 2 = [(1 +/?) + (!- [3/2)m]u 2 



-(l + /3) + (m + l)(2-/3) 



1/2 (9(1 + f3) 2 - 2(m + 4)(2 - /?)(! + 0) + (to + 2) 2 (2 - (3) 2 ) 



-(l + /3) + (m + l)(2-/3) 



2\l/2 



■Jj 



- 1/2 (25(1 + f3) 2 + 2(m - 4)(2 - (3){l + (3) + (m + 2) 2 (2 - (3) 2 ) 



2\l/2 



u 2 . 



(17) 



In the notation of Fliesser et al. pq| , these modes are labeled by three quantum numbers (n,j,m), where m is the 
usual azimuthal quantum number and n and j are two others which distinguish the radial and angular character of 
the modes. The four families we have identified correspond to (0,0, m), (1,0, m), (2,0, to) and (3,0, to), respectively. 
This identification can be checked by considering the limiting values of the analytic result for j3 — > — 1 and (3 — ► 2, as 
presented in Table 1. For example, the 1-modes in Fig. 1(a) tend to zero as (3 — > 2 and to a finite value for (3 — ► — 1, in 
accord with the analytic results. However the finite limiting values of the numerical results deviate from the analytic 
value of y3ro/2, the difference increasing with increasing to (the agreement is exact for the lowest mode in this family 
since it corresponds to the center of mass mode) . The discrepancy between the numerical and analytic results is not 
due to inaccuracies in the numerical calculations, since we have checked that the numerical basis set is sufficiently 
large to yield accurate results for the modes shown. Rather, the differences are real and reflect the fact that our 
calculations are for a finite number of particles, N, while the Thomas-Fermi results correspond to the N — > oo limit . 
Similar deviations appear for the 2-modes in Fig. 1(b). The analytic limits for these modes are -s/3 for (3 — ► 2 and 
y/Zm/2 for (3 — * — 1. The lowest 2-mode in Fig. 1(b), the odd-parity center of mass mode, does extrapolate to \/3 
but the higher modes in this family, to within our numerical accuracy, do not. In other words, the finite size of the 
trapped gas leads to a weak m-dependence of the limiting value, which presumably disappears as N — ► oo. As noted 
by Fliesser et al. |15|, the convergence with TV is more rapid for modes having smaller n values, which is consistent 



with our numerical findings. 

Another feature of the mode spectrum worth noting is a number of instances of anticrossing-type behaviour between 
modes of equal to and the same parity. The most obvious example of this occurs in Fig. 1(a) for the two to = 
modes near [3=1. Similar behaviour is seen for two of the to = odd-parity modes in Fig. 1(b) near f3 — —0.5. 
Of course, there is no avoided crossing of modes with different to values, and there are numerous examples of these 
crossings in Fig. 1. However for real traps which do not possess ideal cylindrical symmetry, one can expect to see 
additional anticrossings associated with the coupling between modes induced by trap imperfections. 



IV. CONCLUSIONS 



In conclusion, we have presented a new and efficient method for solving the Bogoliubov equations for a gas of 
weakly interacting Bose condensed atoms in an axially symmetric magnetic trap. We have used the method to 
calculate the low lying collective mode frequencies and densities over a large portion of the anisotropy parameter 
space. The evolution of the modes as a function of the anisotropy parameter (3 has provided a more complete 
understanding of the relationship of the modes in anisotropic traps to those in the isotropic limit. We have also 
made contact with analytical results obtained in the Thomas-Fermi limit and have identified differences between the 
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finite-TV and infinite- N calculations. We are presently applying our method for anisotropic traps to the problem of 
finite-temperature excitations in the hydrodynamic regime |lq| . 
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Table 1: Limiting values of the mode frequencies, LOi/ui, given in Eq.(|l7j). 
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FIGURE CAPTIONS 

Fig.l: Mode frequencies in units of the average trap frequency u> as a function of the anisotropy parameter (3: (a) 
even- parity modes and (b) odd-parity modes. The symbols code the various m- values: to = (open circle), 
m = 1 (filled square), m = 2 (open triangle), m = 3 (filled circle), m — 4 (open square). The numerical labels 
1-4 identify the different families of modes as discussed in the text. 

Fig. 2: Mode density (left panel) and corresponding equilibrium density (right panel) for values of the anisotropy 
parameter shown. The solid line gives the variation along the x-axis and the chain curve along the z-axis. 

Fig. 3: Contour representation of the mode density of the lowest mode (to = 0, even parity) labeled '3' in Fig. 1(a): 
(a) (3 = —0.5, (b) p = and (c) (3 = 1.4. The shaded region corresponds to negative values of the mode density, 
the unshaded region to positive values. 

Fig. 4: As in Fig. 3, but for the lowest mode (to = 0, odd parity) labeled '4' in Fig. 1(b). 
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